(* Mathematica Package *)

BeginPackage["KramersKronigOscillators`"]
(* Exported symbols added here with SymbolName::usage *)  
ClearAll@ImKramersKronigTransformRe
LorentzianOSC::usage = "LorentzianOSC[\[Omega],{\[Omega]0,A,\[CapitalGamma]}] is the solution to a damped, driven oscillator with analytical real and imaginary parts that are Kramers-Kronig consistent."
VoigtOSC::usage = "VoigtOSC[\[Omega],{\[Omega]0,A,\[CapitalGamma],\[Sigma]}] is the LorentzianOSC convolved with a normal distribution of states"
GaussianOSC::usage = "GaussianOSC[w_, {w0_, A_, s_}] is the KK consistent oscillator with Gaussian imaginary component, this is at present a very slowly evaluating function"
ImKramersKronigTransformRe::usage = "ImKramersKronigTransformRe[ImF, x] calculates the KK transform of the imaginary part of F to get the real part of a KK consistent oscillator.  ImF must be a pure function."

Begin["`Private`"] (* Begin Private Context *)

<< ThinFilmNLSupdater`
CheckUpdate[DirectoryName @ $InputFileName]
 

LorentzianOSC[w_, {w0_ , A_, G_}] := A / ( -(w - w0) - I * G) 


z[w_, {w0_, G_, s_}] := 
	(w - w0 + I * G)/(Sqrt[2] * s);

FadeevaF[z_] := Exp[-z^2] * Erfc[-I * z];

VoigtOSC[w_, {w0_, A_, G_, s_}] := I * Sqrt[\[Pi] / 2] * A / s * FadeevaF[z[w, {w0, G, s}]];


ImKramersKronigTransformRe[ImF_, x_, opts:OptionsPattern[]] := 
	2 / \[Pi] * NIntegrate[xp * ImF[xp] / (xp^2 - x^2), {xp, 0, x, \[Infinity]}, Evaluate[{Method -> "PrincipalValue"} ~Join~ FilterRules[{opts}, Options[NIntegrate]]]];

GaussianOSC[w_, {w0_, A_, s_}, opts:OptionsPattern[]] := I * A Exp[-((w - w0) / s)^2] + ImKramersKronigTransformRe[A (Exp[-((# - w0) / s)^2] - Exp[-((# + w0) / s)^2])&, w, opts];


End[] 

EndPackage[]